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ABSTRACT 

The kinetics of the annihilation process, A + A 0, with ballistic par- 
ticle motion is investigated when the distribution of particle velocities is 
discrete. This discreteness is the source of many intriguing phenomena. 
In the mean field limit, the densities of different velocity species decay in 
time with different power law rates for many initial conditions. For a one- 
dimensional symmetric system containing particles with velocity and ±1, 
there is a particular initial state for which the concentrations of all three 
species as decay as For the case of a fast "impurity" in a symmetric 

background of -|- and — particles, the impurity survival probability decays 

as exp(— const, x In^t). In a symmetric 4- velocity system in which there are 
particles with velocities ±f i and ±f2, there again is a special initial condition 
where the two species decay at the same rate, with a = 0.72. Efficient 
algorithms are introduced to perform the large-scale simulations necessary 
to observe these unusual phenomena clearly. 



1. INTRODUCTION 

In this article, we describe some intriguing aspects of the reaction kinetics in single 
species annihilation, A + A ^ 0, when particles move ballistically with a discrete distri- 
bution of velocities. Unexpected long-time phenomena occur which depend fundamentally 
on the form of the initial velocity distribution. The results discussed here are complemen- 
tary to our earlier work on ballistic annihilation with a continuous distribution of particle 
velocities [1]. For this latter system, the exponents characterizing the decay of the con- 
centration and the typical velocity depend continuously on the form of the initial velocity 
distribution. For discrete velocity distributions, however, the decay kinetics exhibits a fun- 
damentally richer character with fundamental differences in long-time behavior for small 
changes in the initial conditions. 
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Our investigation is also inspired by earlier work by Elskens and Prisch and indepen- 
dently Krug and Spohn who considered the kinetics of the "2-velocity", or "±" model in 
one dimension [2] . Here, the initial velocity distribution of reactants is 

P{v,t = 0) =p+6{v - vo) +p-6{v + vo), 

with p+ + = 1. The spatial distribution of reactants has minimal influence on the 
kinetics as long as the distribution is non-singular. For convenience, we therefore consider 
the distribution to be Poisson in this paper. When p+ > the majority species quickly 
reaches a finite asymptotic limit, while the minority density decays exponentially in time. 
In the interesting situation where the initial densities of the two species are equal, the 
density decays as [2] 

c{t) <x ^/c{0)/vot (1) 

in the long time limit. This relatively slow decay, compared to the rate equation prediction 
of stems from initial density fluctuations. In a region of length L there will typically 
be an imbalance of the order of ^/L in the number of + and — particles. After a time t = L 
has elapsed, only this initial number difference will remain within the region. Therefore 
the local particle number is proportional to y/L, and Eq. (1) follows. Thus the system 
organizes into domains of like velocity particles whose typical size grows linearly in time 
as the reaction proceeds (Fig. 1). 

Consider now a simple and natural generalization to the "3- velocity" model [3] . With- 
out loss of generality, the initial distribution of velocities may be written as 

P{v, t = 0)= p+5{v -v+)+ p-5{v + 1) + poS{v), (2) 

with p^ + p_ + pq — 1. We will primarily focus on the symmetric case where — 1 
and = p^ = P±- The space-time evolution of this system in one dimension for two 
representative values of ip±,Po) is shown in Fig. 1. One of our basic goals is to under- 
stand the time dependence of the mobile and stationary concentrations for different initial 
conditions. Particularly intriguing is the transition from a regime where the stationary 
particles persist, for po > 1/4, to a regime where stationary particles decay more rapidly 
than the mobile particles, for po < 1/4. At a "tricritical" point located at po = 1/4, the 
concentrations of both the mobile and stationary species decay as [4]. While there 

is now a theoretical approach to compute this exponent exactly [5], there is not yet an 
intuitive understanding of this striking behavior. Another intriguing facet of this system 
is the decay of a "fast impurity" , namely, a single particle with velocity -|-1 in system with 
equal concentrations of and —1 particles. By considering the dominant contributions 
to the impurity survival probability, we flnd that this quantity decays asymptotically as 
exp(— const, x In^t). 

Another class of interesting behavior is exemplifled by the symmetric 4- velocity model 
with initial velocity distribution 

p{v, t = 0)= pi{5{v - vi) + 5{v + vi)) + P2i5{v - V2) + 5{v + ^2)), (3) 

with V2 > vi and pi + P2 — 1- According to the rate equations, the more mobile species 
decays as t~'"2/'"i while the less mobile decays as t~^, independent of the initial concen- 
trations. In one dimension, however, either the slower or the faster particles dominate 
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in the long time limit, depending on their relative initial concentrations. At a critical 
value of Pi/p2 which depends on Vi/v2-, numerical simulations indicate that both species 
decay as with a = 0.72. We shall also argue that systems with symmetric velocity 
distributions with n > 4 components exhibit behavior which is characteristic of either the 
3-velocity model, if n is odd, or 4-velocity model, if n is even. Thus we focus primarily on 
the 3- and 4-velocity models as the simplest in the family of symmetric discrete velocity 
models. 

In the next two sections, we discuss the annihilation kinetics of discrete velocity models 
in the mean-field limit. We first treat the conventional rate equations which are based on 
one-dimensional kinematics. The shortcoming inherent in the assumption naturally leads 
us to consider "constant speed" models and the correct d-dimensional kinematics. We find 
rich kinetic behavior which depends on the ratios of initial concentrations, the particle 
radii, and the speeds of the different species. In section 4, we study the kinetics of one- 
dimensional systems. Given the subtle nature of many of our observations, relatively 
efficient and specialized algorithms were developed to provide sufficient data to determine 
the long-time behavior with confidence. For the symmetric 3-velocity model, we investigate 
the decay associated with the tricritical point and the exp(— const, x \v? t) decay of 

the fast impurity problem, in detail. We then consider the kinetics of the symmetric 4- 
velocity model, once again concentrating on the multicritical behavior associated with the 
initial condition where both the fast and slow particles decay at the same rate. 



2. MEAN FIELD THEORY WITH ONE DIMENSIONAL KINEMATICS 

(a) The 3- Velocity Model 

The mean-field rate equations for the 3-velocity model are deceptively simple, but lead 
to relatively complex behavior. For symmetric velocity distribution, the rate equations for 
the concentrations of the left- moving, right-moving, and stationary species, c_(t), c+(t), 
and co(t), respectively, are 

Co = -co(c+ + c_), 

c+ = -c+(co + 2c_), (4) 
c- = -c_(co + 2c+), 

where the overdot denotes time derivative. The numerical factors of 2 reflect the fact that 
the rate of a H — collision is twice that of -f- or — collisions, if we assume that particles 
move only in one dimension. It is in this spirit that the above rate equations are referred to 
as mean-field theory with one-dimensional kinematics. A more complete approach which 
incorporates (i- dimensional kinematics will be outlined in the next section. 

To solve these equations, it is helpful to rewrite the rate equations in terms of ip = 
(c_|_ + C-)/co and = (c_|_ — c_)/co, and the modified time dx = cq dt. This gives, 

^j' + ij = <t^\ 
0' + = 



where the prime denotes differentiation with respect to x. Use of the integrating factors 
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$ = and $ = simplify these equations to 



(6) 



from which it is evident that — $2 _ ^onst. = > 0. Thus the equation of motion for 
W becomes = e~^(^'^ — a^), with solution 



= acoth 



coth 



-1 



*(0) 



- ay 



(7) 



a 



where he have introduced the new time-like variable dy = dx, with y = 1 — e 
monotone increasing function of t. 

To classify the long-time behavior, consider the relative composition triangle p+ + 
P- +po — 1. As indicated in Fig. 2, a given initial condition typically evolves to a "phase" 
where only a single species persists in the long time limit. Consider first the stationary, 
or "0", phase, co(oo) > 0. From the definitions of x and y, the condition co(oo) > 



implies that — e~" ast— >oo, which leads to i/j ~ e~ . Thus the concentrations of 

the mobile species decay exponentially in time in the phase. Interestingly, for p^ = p_ 
stationary particles persist even if the initial concentration of stationary particles becomes 
small. However, the width of this phase becomes vanishingly small in this limit. To 
determine this width, note, from Eq. (4) , that on the boundary where the densities of the 
stationary and positive particles decay at the same rate, the asymptotic solutions of the 
rate equations are c_|_, cq oc 1/t while c_ ~ Thus the boundary between the and + 

phases can be identified by the ratio i/j approaching a finite limit as t — > 00. Since 1 — y 
approaches as t ^ 00, a finite limiting value for i/j = e~^^ requires that the argument 
of the hyperbolic cotangent in Eq. (7) goes to zero. One thereby finds that the width of 
the phase region vanishes as exp(— l/co(0)) as co(0) 0. 

For the symmetric system, a detailed computation gives the asymptotics 



c±(t)~^co(oo)G'(A)e-'^°(°°)*, 



co{t) ~ Co (00) exp G(A)e 



-co(co)t 



(8) 



with G{X) = exp 



dz_^ — z 



V (1 - e"^) ' ^ = 2c±(0)/co(0), and with the final den- 



sity of stationary particles given by 

co(oo) = co(0) e-^'^iW/^oW^ 



(9) 



Thus while a residue of stationary particles always persists, their concentration is astro- 
nomically small if the initial concentration of stationaries is small. 

For the asymmetric 3-velocity model, the corresponding rate equations are 



Co = -co{vc+ + C-), 

c+ = -c+{vco + (1 + v)c-), 

c_ = -c_(co + (1 + v)c+). 



(10) 
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While we are unable to solve these equations, we can readily find the asymptotic behavior. 

In the phase, the rate equations for the mobile particles have the asymptotic solutions 
c_|_(t) ~ g-ti'co(oo) g^j^^ C-{t) ~ g-tco(oo)_ Similar exponential decays characterize the 
behavior of the minority species in the other two phases. On the separatrices, however, 
two species decay at the same rate while the minority species decays faster. For example, 
on the separatrix between the and + phases, one has cq ~ c+ ^ c_. Substituting this 
into Eq. (10) yields the asymptotic solution c_|_,co ^ 1/vt while c_ ~ 

t-i-2/v_ Similarly, 

on the separatrix between the and — phases c_, Cq — 1/t and C-|- ~ 

The long-time persistence of stationary particles also occurs in a general 2m + 1- 
component model with velocities = vq < vi < . . . < Vm and corresponding concentrations 
co(t), ci(t), . . . , Cm{t). The rate equations for these concentrations are 

m 

Co = -2co'^VjCj, 

(11) 



k — 1 m 



j=l j=k 

Introducing x = j^dt' Ci{t') and the dimensionless concentrations 0/- = C}./cq^ we obtain 
a closed system of equations for (j)k{x) and an additional equation for cq{x) 



d\n(f) 



dx 
din Co 



dx 



k-i 

= -Vk -'2^{vk - Vj)(j)j, k = l,...,m 

(12) 



Since Vk > Vj for k > j, the first of Eqs. (12) gives (f)k{x) < 0A;(O)e Substituting these 
into the equation for cq yields 

( ) ^ 

^^Tw^^-^T.Mm-e-'^n (13) 

which immediately leads to the lower bound for final density of stationary particles 



co(oo) > co(0) exp 



co(0) ^ 



E^^(o) 



(14) 



For the 3- velocity case this bound is exact and coincides with Eq. (9). Thus stationary 
particles always survive in the symmetric (2m + l)-velocity model, although the final 
residue is vanishingly small when their initial concentration is small. 
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(b) The 4- Velocity Model 

For the symmetric 4-velocity model, denote by ci and C2 the concentrations of species 
with velocities ±vi and ±V2, respectively. Without loss of generality, let V2 > vi and set 
vi — 1,V2 — V. The rate equations for this system are, 

Cl = —2c? — 2vciC2, 

C2 = —2vc2 — 2vCiC2, 

with asymptotic solution 

ci{t)^t-\ C2{t)^t-\ (16) 

Thus the faster species decays non- universally as This asymptotic behavior is reached 
only at very long times, however, when the initial velocities are nearly identical. Essentially 
the same equations were solved in the context of heterogeneous diffusive single-species 
annihilation in which particles have different diffusion coefficients [6] . 

Analogous behavior occurs in the symmetric 2m-velocity model with concentrations 
ci(t), . . . , Cm(t) and speeds ±Vj, with vi < . . . < Vm- The rate equations are 

(k—l m \ 

Vk^Cj + ^VjCj . (17) 
3=1 3=k J 

Introducing now x = 2 dt' ci(t') and (j)k = c^/ci, = 1, we obtain 



—^^^-'^i'"k-Vj)(l)j, 2<k<m, 



d In Cl 

= -Z^V3<l>j- 



(18) 



dx 

3 = ^ 



As in the (2m-|-l)-velocity model, one can straightforwardly derive (pki^) < (j)k{0)e~^^''~^^^^ • 
This, together with the equation for ci and the relation 0i = 1, show that ci ~ e~^^^. 
Combining this result with the definition of x proves that a; — > oo as t — > oo. It therefore 
follows that in the long time limit (j)kix) ~ e"^'''^"'"^)^. Re-expressing this in terms of Cj{t) 
leads to 

ci-t-\ C2~t-^^/^S c^~t-^-/^i. (19) 

Thus a more mobile species k decays non-universally with an associated exponent equal 
to the velocity ratio Vk/vi. 

3. MEAN FIELD THEORY WITH D-DIMENSIONAL KINEMATICS 



We now generalize the rate equations to account for d-dimensional kinematics. This 
should be viewed as the "true" mean-field theory of ballistic annihilation. It is convenient 
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to consider first the kinetics of a single impurity in a background of particles moving at 
the same speed. From this, the general mean-field theory follows naturally. 

(a) The Impurity Limit 

Let background particles of radii R move with identical velocities v which are uni- 
formly distributed in angle, i. e., P(v,t = 0) = 5(|v| — vq). If we temporarily neglect the 
annihilation events among background particles, then by an elementary mean-free path 
argument, the (infinitesimal) concentration of the impurity species with radius Rj varies 
as 

ci = -cicvnd-i{R + Ri), (20) 

where Qd-i{R) is the volume of a sphere of radius R in d — 1 dimensions. 

If the impurity moves with velocity w, then the decay rate must be averaged over all 
directions of relative velocities, v — w. This leads to 

, C de isinef-^^l + e2 + 2ecos^ 

Jq dO (sm^)" ^ (21) 

= -cicxvVLd-i{R + Ri) X T{e). 

with e = w/v. To find the impurity concentration c/ we must first determine the back- 
ground concentration c. Since a background particle can be considered as an impurity of 
radius R moving with velocity v, we apply Eq. (21) to obtain c = — v f2d_i(2i?) J^(l), 
with asymptotic solution c{t) ~ [vf2d-i(2-R).F(l)t]~-'^. Using this in Eq. (21) gives 

crr^t , with a = a{R„e)= ■ (22) 

For example, when there is a stationary particle in a uniform background of moving par- 
ticles, c/ ~ t-«o(d)^ with 

a„(rf)=^(0)mi) = A|Li_2l, (23) 

and F is the gamma function. Note ao{d) is rational for odd dimensions and transcendental 
for even dimensions: ao{l) = 1, ao{2) = tt/A, ao{3) = 3/4, etc.. Interestingly, q;o(oo) = 
l/-\/2, which can be understood by noting that in the limit d — > oo, two arbitrary particles 
always move orthogonally with relative velocity vy/2. The expression for q;o(oo) (Eq. (22)) 
essentially involves the inverse of this factor. 

(b) Stationary and Moving Species 

Consider now the (i-dimensional analog of the symmetric 3-velocity model in which 
there is a finite initial concentration of stationary particles of radii Rq and mobile particles 
of radii R all moving with speeds v. The appropriate rate equations for the corresponding 
concentrations co{t) and c{t) are similar to those for the above impurity limit except that 
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one must account for the influence of stationary particles on moving particles. The rate 
equations become (compare with Eqs. (4) with C-(- = c_) 



Co = -vQd-iiR + Ro)coc, 
c = -vQd-iiR + Ro)[coc + XdC 



(24) 



where = J-'{l)Qd-i{^R)/^d-iiR + Ro)- Introducing the modified time variable, T = 

gives the linear equations 

dco 

^ = -co - A,c, 

which can be readily solved to give 

co(T) = co(0)e-^, 

e-^-e-^<^^ (26) 
c(r) = c(0)e-^'^^-co(0) " ^^^^ ■ 

For Ad > 1, the concentration of the mobile species decays exponentially in real time 
t while stationary particles always persist: 



co(oo) = Co(0) 



1 + (Ad - 1) 



c(0) 
co(0) 



(27) 



For Ad < 1, there are three different regimes. For small initial concentration of 
mobile species, c(0) < 1/(2 — Ad), the density of moving particles decays exponentially in 
time, while stationaries persist with a residue still given by Eq. (27). At the critical point, 
c(0) = 1/(2— Ad) andco(O) = (l — Ad)/(2 — Ad), both species decay as Finally, for c(0) > 
1/(2 — Ad), both species decay as distinct power laws in time: c{t) ~ 3> co(t) ~ 
Amusingly, these same three qualitative regimes occur in the one-dimensional 3-velocity 
system. For example, in three dimensions, A3 = \/24( a+^q )^i hence the threshold between 
different behaviors, A3 = 1, occurs when Rq/R = 3.42673. Thus the stationary particles 
always persist if their relative size is small, Rq/R < 3.42673, while for Rq/R > 3.42673, 
stationary particles may disappear if their initial concentration is small enough. 

(c) The 2-Speed Model 

For the d-dimensional 2-speed model, with Cj,Rj, and Vj the concentration, radius, 
and speed of j^^ species, j = 1,2, with e = Vi/v2 < 1, the rate equations are 

6i = -cf-MC,c. ^^^^ 

C2 = -/XC1C2 - 1^C2, 

where we have set vif2d_i(2i?i).F(l) = 1 by rescaling the time and have introduced 

1 Qd_,iR, + R2) Tje) 1 »d-i(i^2) 
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Eq. (28) give rise to three asymptotic behaviors which depend on the relative magni- 
tudes of fi and v: 



Ci{t)^t-\ C2(t)~t"'', 1<U<H, 

ci{t)^t-^/\ C2{t)^t-\ z/ <//<!, (30) 

Ci{t)^t~'^, C2{t)^t~^, V<IJL<l,IJi<l<V. 

In the two remaining cases of 1 < i/ < and v <1 < ji the asymptotic behavior depends 
also on initial concentrations. The first and the second asymptotics are realized when 
^ c^(o) ; respectively, while the third asymptotics occurs if = ^^{sj' Parentheti- 
cally, for an arbitrary number of mobile species of equal radii in three dimensions, a similar 
analysis gives 

c,W~r-, .. = ^ + g, (31) 

where vi is the smallest velocity (therefore, jii = 1). Thus the least mobile species decays 
as t~^, while the more mobile species decay non-universally. 

4. KINETICS IN ONE DIMENSION 

(a) Geometric Approach for the Symmetric 2- Velocity Model 

For completeness and to provide a framework to discuss the 3- and 4-velocity models, 
we first give a geometric derivation for the decay of the concentration in the 2-velocity or 
± model when the initial concentrations of the two species are equal. This approach is 
based on the equivalence between the kinetics of the particle system and the smoothing of 
one-dimensional stepped interface (Fig. 3). In this mapping, a right- moving particle in the 
± model, is equivalent to an "up" step in the corresponding interface. This up step moves 
to the right at the same speed as the initial particle. Similarly, a left-moving particle is 
equivalent to a left-moving "down" interface step. An annihilation event in the particle 
system corresponds to the disappearance of a tier in the interface. Clearly, the collision 
partner of a given up step is the first down step to the right which is at the same height 
as the initial step. 

Since up and down steps steps are uncorrelated and occur with equal probability, the 
probability that the initial up step (defined to be particle 0) collides with its (2n -|- 1)^* 
neighbor is exactly equal to the first passage probability for a random walk to return to the 
origin at 2?i + 2 steps. Thus the initial right-moving particle annihilates with a left-moving 
(2n -|- 1)^* particle with probability [7] 

/n = 2-^"-^^/^ (32) 
n\[n -\- 1)! 

Asymptotically, this collision probability varies as n""^/^ as n — > oo. Consequently, the 
probability that the initial particle survives potential encounters up to the 2n*'^ neighbor, 
ySji = 1 — Ylm'Knfn'i leads to a survival probability which decays as Sn ~ n"-*^/^, or 
c{t) ~ t~^^'^ in the continuum limit. 
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(b) The 3- Velocity Model 



For the 3- velocity model, there does not appear to be a similar geometric construction 
to help determine the kinetics. We therefore resort to qualitative arguments, as well as 
numerical simulations, to determine the long-time behavior. In one dimension, we expect 
the kinetics to be different from mean-field predictions because of the tendency of like 
velocity particles to cluster, as observed in space-time graphs (Fig. 1). For concreteness, 
we focus on the symmetric system where pj^ — p- = p± and first investigate whether 
stationary particles persist for any value of po by numerical simulations. Because of subtle 
crossover effects, a direct molecular dynamics approach is inadequate to yield accurate 
results and we therefore developed a more efficient approach in which all collision partners 
and corresponding collision times are identified at the outset. 

In this algorithm, stationary or right-moving particles are placed on a stack (first in, 
last out) as they are initially created. When a left-moving particle is created, its collision 
partner is determined immediately, since this partner is necessarily one of the particles from 
the already existing stack. (There is a particular case in which a negative velocity particle 
is deposited when the stack is empty. This exits the system, since free boundary conditions 
are employed. To ensure that this effect does not give spurious results, only particles from 
the middle half of the system are considered) . The determination of the collision partner of 
the left-moving particle is accomplished by straightforward comparisons. If the uppermost 
particle on the stack moves to the right, then it is the collision partner, and the collision 
time is recorded. On the other hand, if the "last" particle on the stack has zero velocity, one 
must compare the collision time between this last particle on the stack and the left-moving 
particle, and the collision time between the last particle with "earlier" outgoing particles 
from the stack. All collision partners and corresponding collision times are determined, 
up to and including the collision time of the initial left-moving particle. Right-moving or 
stationary particles for which the collision time is determined are then removed from the 
stack. Prom the stored array of collision times one determines c±{t) and co(t) by counting 
the number of particles of a given species that survive at that time. 

With this method we simulated 5 x 10^ particles to 10^ time steps in approximately 
30 cpu seconds on a DEC/AXP 3000/400 workstation. Our numerical results are typ- 
ically based on 100 realizations at each initial concentration. These simulations reveal 
the following basic results (Fig. 4): For pq < 1/4, co{t) ~ 1/t and c±{t) ~ The 
crossover to the asymptotic behavior becomes progressively more gradual as po — > 1/4 
from below and there is a substantial time range for which c±{t) and co{t) decay at nearly 
the same rate before the final asymptotics is reached. (This is the primary reason for 
the erroneously-reported non-universal behavior based on data from direct and much less 
extensive molecular dynamics simulations [1].) Exactly at po = 1/4, the data indicate 
that both c±{t) and co{t) decay as t~'^, where extrapolations of local slopes of neighbor- 
ing points in figure 4(b) give a = 0.665 (Fig. 5). When the time difference between these 
points is the factor 1.2^^ ^ 46, a smooth sequence of local exponents is obtained. However, 
similar results are obtained for different choices of time delay factors. The trend in the 
local exponents suggests, in fact, that ct = 2/3, a result which has now been obtained by an 
exact calculation [5]. Additionally, we estimate the relative amplitude c±{t)/co{t) = 1.17. 
For Po > 1/4, co{t) saturates to a finite limiting value which appears to be proportional to 
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(po — while c±{t) decays faster than a power law in time. Based on these results, the 
phase diagram shown in Fig. 6 is inferred. 

The location of the tricritical point, where all three species decay at the same rate, 
may be found by the following heuristic argument [8]. Since half the stationary particles 
react with + particles, the fraction of + particles available to react with — particles is 
simply — ^pq. This is proportional to the number of H — annihilation events per unit 

length, Similarly, the relative number Pq- of 0— annihilation events per unit length 

is equal to |po- Now we assume that /Pq- = 2, based on the expectation that 

the relative number of annihilation events is proportional to the relative velocities of the 
collision partners. Combining the resulting relation, — |po = POj with the normalization 
condition, 2p+ + pq = 1, we find the location of the tricritical point to be po = 1/4, p+ = 
p_ = 3/8. 

It is straightforward to generalize this argument to the asymmetric velocity distribu- 
tion ,0, —1). Since the + — symmetry is now broken, we must consider separately 
the ratios, P+_/Po- and P_+/Po+- Following the same considerations as in the symmetric 
case, these reaction numbers are 

P+- ^ P+ - vpo/i^ + v) 
Po- ~ Po/{l + v) 

P-+ ^ P- -Po/i'^ + v) 
Po+ vpo/{l + v) 

Together with the normalization condition, we find, for the initial concentrations at the 
tricritical point, 

f» = j. f+ = j(i + Ti^). f- = K'^TT^)- '''' 

While Eq. (34) involves uncontrolled approximations, especially since symmetry con- 
siderations no longer apply, it is relatively accurate. For example, simulations with v = 2 
suggest that the tricritical point is located at po ~ 1/4 and p+ ~ 0.402, compared to 
= 5/12 = 0.416 from Eq. (34). Generally when v varies from to oo, the tricritical 
point moves along the straight line parallel to the base of the composition triangle. At 
the extreme limits oi v = and v = oo the 3-velocity model degenerates to the 2-velocity 
model for which we know that Eq. (34) is exact. The available numerical evidence also 
supports the apparent general equality po = 1/4. Another interesting question in the gen- 
eral asymmetric case is whether the three coexistence curves really coalesce at one point 
or whether there is an open region where all concentrations decay at a similar rate. Nu- 
merics cannot answer such a question definitively, but the evidence appears to favor the 
hypothesis of one single tricritical point where all three coexistence lines merge. 

Numerical results also suggest that near the tricritical point of the symmetric model, 
Po = l/4,p± = 3/8, the long time kinetics depends on a single scaling variable ^ = tA^ 
where A = ~ 1/4 is assumed to be small. The scaling assumption for the concentrations 
gives 

co{t, A) ~ t-^/^Co{tA^), c±{t, A) ~ t-^/^C±{tA^), (35) 



= l + v, 
1 + v 



(33) 
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with Co(0 ^±(0 finite at ^ = to reproduce the decay at the tricritical point. 

This impUes that Co(t, A) ~ A^, as is observed numerically. For A < the scaling 
predictions agree with simulations if Co(0 ~ ("0""*^^^ ^^d C±{^) ~ (—0"*^^^ as ^ — —oo. 
Thus scaling provides the A dependence of the asymptotic behavior for A < 0, namely, 
co(^) ~ (^A)~^ and c±(t) ~ (A/t)^/^. Additionally, the crossover time from tricritical 
behavior to the final asymptotics scales as A~^. 

Since the concentration of stationary particles saturates to a finite value for A > 0, 
Co(0 ~ ^^^^ as ^ ^ oo. In this regime, c±{t) will eventually decay rapidly, presumably 
exponentially, leading to a similar decay of C±(^) as ^ — > oo. This has a remarkable 
consequence for the spatial distribution of immobile particles. Since the crossover time is 
of order A~"^, c±{t) should decay asymptotically as exp(— tA^). The residual concentration 
of the immobile particles, however, goes as A^. The two facts appear contradictory, since 
the decay law for c±{t) implies that immobile-free intervals of length of order A~^ must be 
reasonably frequent. The resolution of this difficulty is presumably that immobile particles 
are distributed on a fractal set on length scales up to A~^. A consistent value for the fractal 
dimension of this set is 1/3, since this value indeed implies that there are A"-*^ particles in 
an interval of length A~^. This observation has been verified numerically. 

(c) The Impurity Limit of the 3- Velocity Model 

When the concentration of one species is negligible while the other two species have 
equal concentrations, it is possible to determine the asymptotic survival probability of the 
"impurity". These results are worth emphasizing, both for methodological interest and 
because of their unusual nature. First consider a stationary impurity in a background of ± 
particles. By the mutual annihilation of ± particles, their density decays as c±{t) ~ 
On the other hand, a stationary particle survives only if it is not annihilated by particles 
incident from either direction. Since each of these two events is independent, it follows that 
co{t) ~ c±{t)'^ ^ in the limit po <ti p±. This same argument continues to apply if the 
impurities are "slow", i. e., their velocity w satisfies \w\ < 1. The full survival probability 
is again a product of single-sided survival probabilities, each of which decays as t"-*^/^. 

More precisely, svippose that the impurity starts at the origin with velocity w and 
that there is a Poisson distribution of initial separations of the background particles. The 
impurity survival probability at time t, S{t, w), equals the product of left- and right-sided 
survival probabilities, S-{t,w)S+{t,w). Since S-{t,w) = S+{t,—w), we only need to 
compute S+{t,w). This one-sided survival probability is given by 

oo 

S+it,w) = J2fnP[^2n+l > il + w)t]. (36) 
n=0 

This is simply the sum over all collisions partners of the probability that the collision time 
between the impurity and each of its potential collision partners is greater than t. For a 
Poisson distribution of interparticle separations with unit density, the probability that a 
left-moving (2n -|- 1)*^* particle is located at X2n+i > (1 + '>^)'t is 

P[x2n+i > (1 + w)t] = / -—e-'dz. (37) 

J{i+w)t (2n)! 
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Combining Eqs. (32), (36), and (37) yields [9] 



S+it,w) = l- [ 



il+w)t 



dz 




z 



(38) 



e-^'+^^\lo[il + w)t]+h[il + w)t]). 



where Ii is the modified Bessel function. Thus the complete survival probability is 



in agreement with the rough argument given above. 

Particularly intriguing behavior occurs in the complementary "fast" impurity limit 
with a vanishingly small initial concentration of + impurities in a background of equal 
concentrations of — and particles. By a Galilean transformation, this system can be 
viewed as an impurity of velocity +3 in a H — background. More generally, we consider 
the decay of a fast impurity with speed |w| > 1 in a H — background. An asymptotic 
argument suggests that this survival probability decays slower than an exponential but 
faster than a power law in time. 

The basis of this argument is to consider a subset of configurations which gives the 
dominant contribution to the impurity survival probability, but which are sufficiently sim- 
ple to evaluate. For the impurity to survive to time t, the background ± particles must an- 
nihilate among themselves up to this time. On a space-time diagram, these self-annihilation 
events appear as a sequence of isosceles triangles which do not extend to the world line of 
the impurity. We posit that the dominant contribution to the survival probability stems 
from a sequence of systematically larger self-annihilation triangles which just miss the im- 
purity world line (Fig. 7). From basic geometry, the base of the n*^ triangle, Xn, equals 
xqP'^, with P = {w + l)/{w — 1). Here the separation between successive triangles is ne- 
glected, an approximation which is valid as the number (and size) of the triangles becomes 
large. Under this assumption of abutting triangles, the number of such triangles that 
comprise the self-annihilation sequence up to time t is ~ ln{t/xo)/ ln(3. 

By construction, collision partners in the background define the sides of an individ- 
ual self-annihilation triangle which encloses more local self-annihilation events. If these 
collision partners arc separated by a distance x„, the existence probability of the self- 
annihilation triangle is (x^/xq)"'^/^ as n oo. The impurity survival probability S{t) 
is therefore the product of occurrence probabilities for the sequence of self-annihilation 
triangles up to time t, leading to 

N 



S{t, w) = e-^\lo[{l + w)t] + /i[(l + w)t]) X (/o[(l - w)t] + /i[(l - w)t]), 



2 



(39) 



t 



-1 



7r\/l — w'^ 




(40) 



n 



exp(-ln^(t/a;o)/|ln/5). 
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This result should be accurate when the number of self-annihilation triangles is large, a 
situation which occurs when w > 1. In this limit, it is impractical to directly simulate 
the large systems needed to confirm the above prediction. We have therefore written a 
program which tracks only the impurity and the potential collision partners at any given 
stage of the reaction. Particles which are part of local self-annihilation events need not, 
and are not considered, leading to negligible CPU requirements. With this method we 
have simulated 10^ realizations of an impurity with w as small as 1.0002. In this case, 
the mean impurity lifetime is approximately 340, but there are a few configurations with 
substantially longer lifetimes. For w < 1.005, a plot of the logarithm of the survival 
probability versus z = In^ t/ln/S exhibits good data collapse and linear behavior over a 
substantial range (Fig. 8). 

From S{t) given in Eq. (40), the mean impurity lifetime, (t) = J°° S{t') dt', is readily 
computed to be proportional to {■^jzij)'^^^- While this is in excellent agreement with numer- 
ical results, the use of the asymptotic form for S{t) even for t of order unity in the integral 
for {t) has yet to be justified. One additional amusing feature is that for — > 1+, the prob- 
ability that the impurity annihilates with a right-moving particle vanishes as (ty — 1)^/"^. 
When w = 1, the impurity becomes one of the right-moving background particles which 
can only annihilate with left-moving particles. 

(d) The 4- Velocity Model 

We focus on the symmetric system with particles of velocities ±Vi and ±V2 and relative 
concentrations pi/2 and P2/2, respectively. According to the mean-field description, the 
more rapid particles disappear more quickly for any pi,p2 > 0. However, we find three 
regimes of behavior in one dimension (Fig. 9). Depending on e = ^1/^2, there is a critical 
value Pi(e) such that for pi < Pi(e), ci{t) ~ while C2{t) ~ t~^^'^. The crossover to 
these asymptotic behaviors sets at progressively later times as pi approaches Pi{e) from 
below. The converse behavior occurs for pi > Pi(e). Roughly speaking, in these two cases 
the system reduces to the 2-velocity model as the minority species disappears. However, 
when pi = Pi(e) both species decay at the same rate. Based on extensive simulations, 
this decay appears to be a power law, with a = 0.72 ± 0.01 (Fig. 9). This value is 
obtained by performing a least-squares fit to the data on a double logarithmic scale in the 
time range where linear behavior is most evident, typically for 10^ < t < 10^. The error 
estimate is based on the variation of the exponent values for systems with pi within 0.01 
of the critical value. 

For detecting this power law behavior, direct molecular dynamics is once again inade- 
quate. We therefore developed a more efficient algorithm which can be viewed metaphor- 
ically as "pick up sticks" . This approach is applicable for any initial velocity distribution 
with non-degenerate collision times. We first identify the collision times of all nearest- 
neighbor pairs and then sort them in ascending order by a standard 0{N In N) algorithm 
[10]. Next, these near-neighbor collision times are sequentially added to the list of true 
collision times if a consistency criterion, to be specified below, is satisfied. At each storage 
event, the particle pair (n, n + 1) associated with the underlying collision is removed from 
the system. Correspondingly, the collision times associated with (n — 1, n) and {n + l,n + 2) 
must be discarded, while the collision time associated with the new nearest-neighbor pair 
(n — 1, n-|-2) is computed. If any of these three collision times is smaller than the next near- 
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neighbor collision time in the previously sorted list, it is necessary to re-sort the current 
list of collision times before continuing with the sequential storages of true collision times. 
Since the largest nearest-neighbor collision times will never be reached before re-sorting is 
necessary, the sorting is performed only on a small fraction of the smallest of these collision 
times at each stage. With this method we simulated 2.5 x 10^ particles to 10^ time steps 
in of the order of 70 cpu seconds on a DEC/AXP 3000/400 workstation. Our results are 
typically based on 50 realizations at each initial concentration. 

The location of the critical point for the symmetric 4-velocity model may also be 
estimated by the same approach that was applied for the 3-velocity model. Let Pjj be 
the number of annihilation events between -|- and — particles of type j, j = 1,2, and let 
Pijj ^ 7^ be the number of annihilation events between say a -|- particle of type i with 
both + and — particles of type j. There are two "mass" conservation laws, P22 + -P21 p2 
and Pii + P12 cc pi, as well as the symmetry condition P21 = Pi2- We further assume 
that P22/P21 — 2, indicative of the fact that the relative velocity between -|- and — V2- 
particles is equal to 2v2 while the average relative velocity between say a -|-V2-particle and 
an arbitrary ±fi-particle is equal to V2- Analogously, one also has P11/P12 = 2v\/v2- 
Combining these relations, the ratio of initial densities at the critical point is 



pi 1 -I- 2i;i/i;2 

Eq. (41) reproduces the expected results in two extreme limits of vi — V2, where it gives 
P2 — Pi, while for vi — 0, the 3-velocity tricritical point is reproduced. For two particular 
cases that were simulated extensively, Eq. (41) gives P2/P1 = 1-5 for e = 2 and P2/P1 = 
1.065 for e = 1.1, while the corresponding estimates from simulations are 1.6 and 1.20. Thus 
while Eq. (41) reproduces the correct qualitative trend for the location of the tricritical 
point, it is less accurate than the corresponding prediction of the 3-velocity model. 

5. SUMMARY 



Ballistic annihilation when the particle velocities are drawn from distributions is a 
problem that appears ripe for further exploration. There is a rich array of phenomena 
in which the underlying discreteness of the particle velocities is crucial. We have focused 
on simple and generic situations to elucidate these features. The mean- field description of 
ballistic annihilation is deceptively simple but leads to rather complex behavior. Particular 
attention was paid to distinguishing between the case of one-dimensional kinetics (which is 
the basis of conventional mean-field theory) and (i-dimensional kinetics, in which averages 
over all particle directions in d dimensions is accounted for. Since the "number" of species is 
infinite, it is difiicult to imagine that large-scale single-species heterogeneities could form. 
The absence of such spatial organization suggests that the mean-field approximation is 
applicable when d > 1. 

In one dimension, the relative initial abundance of various species fundamentally de- 
termines which species predominates in the long-time limit. For the symmetric 3-velocity 
model there is tricritical behavior in which all three species decay as t~^/^ when p± = 3/8 
and po = 1/4. In the fast impurity limit, the survival probability of the impurity decays as 
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exp(— const, x In^t). For the 4-velocity model, mean-field theory predicts that the faster 
species decays with a non-universal power law, while the slower species always decays as 
1/t. However, in one dimension, it is again the relative initial concentrations which deter- 
mine whether the more rapid or the slower species can dominate asymptotically. At the 
threshold between these two regimes, the densities of all species appears to decay at t~'^, 
with a ^ 0.72. 

It is important to mention that Piasecki and Droz et al. have very recently developed 
a powerful analytic method to solve for the kinetics of one-dimensional of ballistic anni- 
hilation models exactly. In particular, they obtain the decay of the density at the 
tricritical point of the symmetric 3- velocity model that we inferred from simulations. Their 
method also appears to be applicable to general velocity distributions. However, even in 
the three velocity model, the construction of explicit solutions from their formalism is a 
formidable task. Thus it still would be desirable to develop either continuum approaches 
or other analytical methods that would provide better intuitive insights into the intriguing 
qualitative features of ballistic annihilation. 
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FIGURE CAPTIONS 



Figure 1. Space-time representation of particle trajectories in the symmetric 2- and 3- velocity 
models, where particles move with velocity or ±1. Shown are: (a) po — 0.0, p± = 
0.50 and (b) po = 0.25, p± = 0.375. 

Figure 2. "Phase" diagram of the symmetric 3- velocity model in mean- field limit within the 
relative composition triangle defined by triangular region p+ + p- + po = 1. The 
regions marked by -|-, — , and are phases where only positive, negative, or stationary 
particles, respectively, persist in the long time limit. Along the boundaries between 
-|-0 or —0, the concentrations of the two competing species decay as t~^, while the 
minority species decays as t~^. The width of the phase region is vanishingly small 
as Po — > 0. 

Figure 3. Equivalence between the time evolution of the 2- velocity model and a one dimensional 
"wedding cake" interface. The collision partner of the initial particle is indicated in 
both the particle and interface representations by the small arrows. 

Figure 4. Representative simulation results for the symmetric 3-velocity in one dimension. Shown 
are double logarithmic plots of c+{t) {+), c_(t) (V), and co(t) (o) versus t for: (a) 
Po = 0.10, p± = 0.45, and (b) at the tricritical point po = 0.25, p± = 0.375. Each 
successive data point represents in increase in t by a factor of 1.2. 

Figure 5. Local exponents for C{t) = C-^-{t) + C-{t) (+) and Co{t) (circles) versus 1/t at the 
tricritical point. These exponents are based on the slopes of points separated by 
T = 1.2^-'^ « 46 from the data of figure 4(b). 

Figure 6. "Phase" diagram of the symmetric 3-velocity model in one dimension. Along the 
dashed line, c±{t) ~ while co{t) ~ t~^. At the tricritical point (circle), all 

species decay as Along the solid lines, the nature of the decay is unknown, 

except very close to the extrema which correspond to the "fast impurity" problem 
(see text). 

Figure 7. World line of a fast impurity in a background of equal concentrations of background ± 
particles. Successive triangles of self- annihilating background particles are indicated. 

Figure 8. Plot of the survival probability of a fast impurity. Shown is In S{t) versus In^ t/\n/3 
for - 1 = 0.0005, 0.001, and 0.002 (+, A, and x, respectively). 

Figure 9. Representative simulation results for the symmetric 4-vclocity in one dimension. Shown 
are a series of double logarithmic plots of ci{t) (□) and C2{t) (o) versus t for: (a) - (c) 
vi = 1 and V2 — 1.1 for pi = 0.40, 0.45, and 0.50, and (d) - (f) vi — 1 and V2 — 2 for 
Pi = 0.35, 0.38, and 0.40. Each successive data point represents in increase in t by a 
factor of 1.5. 
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